Lyapunov Function Search

Adapted from: SOSTOOLS' SOSDEMO2 (See Section 4.2 of SOSTOOLS User's Manual)

using DynamicPolynomials
@polyvar x[1:3]

(DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[x₁, x₂, x₃],)

We define below the vector field $\text{d}x/\text{d}t = f$

f = [-x[1]^3 - x[1] * x[3]^2,
     -x[2] - x[1]^2 * x[2],
     -x[3] - 3x[3] / (x[3]^2 + 1) + 3x[1]^2 * x[3]]

3-element Vector{MultivariatePolynomials.RationalPoly{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}, DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}}}:
 (-x₁x₃² - x₁³) / (1)
 (-x₂ - x₁²x₂) / (1)
 (-4x₃ - x₃³ + 3x₁²x₃ + 3x₁²x₃³) / (1 + x₃²)

We need to pick an SDP solver, see here for a list of the available choices. We use SOSModel instead of Model to be able to use the >= syntax for Sum-of-Squares constraints.

using SumOfSquares
using CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver);

We are searching for a Lyapunov function $V(x)$ with monomials $x_1^2$, $x_2^2$ and $x_3^2$. We first define the monomials to be used for the Lyapunov function:

monos = x.^2

3-element Vector{DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}:
 x₁²
 x₂²
 x₃²

We now define the Lyapunov function as a polynomial decision variable with these monomials:

@variable(model, V, Poly(monos))

\[ ({\_}_{1})x_{3}^{2} + ({\_}_{2})x_{2}^{2} + ({\_}_{3})x_{1}^{2} \]

We need to make sure that the Lyapunov function is strictly positive. We can do this with a constraint $V(x) \ge \epsilon (x_1^2 + x_2^2 + x_3^2)$, let's pick $\epsilon = 1$:

@constraint(model, V >= sum(x.^2))

\[ ({\_}_{1} - 1)x_{3}^{2} + ({\_}_{2} - 1)x_{2}^{2} + ({\_}_{3} - 1)x_{1}^{2} \text{ is SOS} \]

We now compute $\text{d}V/\text{d}x \cdot f$.

using LinearAlgebra # Needed for `dot`
dVdt = dot(differentiate(V, x), f)

\[ \frac{(-8 {\_}_{1})x_{3}^{2} + (-2 {\_}_{2})x_{2}^{2} + (-2 {\_}_{1})x_{3}^{4} + (-2 {\_}_{2})x_{2}^{2}x_{3}^{2} + (-2 {\_}_{3} + 6 {\_}_{1})x_{1}^{2}x_{3}^{2} + (-2 {\_}_{2})x_{1}^{2}x_{2}^{2} + (-2 {\_}_{3})x_{1}^{4} + (-2 {\_}_{3} + 6 {\_}_{1})x_{1}^{2}x_{3}^{4} + (-2 {\_}_{2})x_{1}^{2}x_{2}^{2}x_{3}^{2} + (-2 {\_}_{3})x_{1}^{4}x_{3}^{2}}{1 + x_{3}^{2}} \]

The denominator is $x[3]^2 + 1$ is strictly positive so the sign of dVdt is the same as the sign of its numerator.

P = dVdt.num

\[ (-8 {\_}_{1})x_{3}^{2} + (-2 {\_}_{2})x_{2}^{2} + (-2 {\_}_{1})x_{3}^{4} + (-2 {\_}_{2})x_{2}^{2}x_{3}^{2} + (-2 {\_}_{3} + 6 {\_}_{1})x_{1}^{2}x_{3}^{2} + (-2 {\_}_{2})x_{1}^{2}x_{2}^{2} + (-2 {\_}_{3})x_{1}^{4} + (-2 {\_}_{3} + 6 {\_}_{1})x_{1}^{2}x_{3}^{4} + (-2 {\_}_{2})x_{1}^{2}x_{2}^{2}x_{3}^{2} + (-2 {\_}_{3})x_{1}^{4}x_{3}^{2} \]

Hence, we constrain this numerator to be nonnegative:

@constraint(model, P <= 0)

\[ (8 {\_}_{1})x_{3}^{2} + (2 {\_}_{2})x_{2}^{2} + (2 {\_}_{1})x_{3}^{4} + (2 {\_}_{2})x_{2}^{2}x_{3}^{2} + (-6 {\_}_{1} + 2 {\_}_{3})x_{1}^{2}x_{3}^{2} + (2 {\_}_{2})x_{1}^{2}x_{2}^{2} + (2 {\_}_{3})x_{1}^{4} + (-6 {\_}_{1} + 2 {\_}_{3})x_{1}^{2}x_{3}^{4} + (2 {\_}_{2})x_{1}^{2}x_{2}^{2}x_{3}^{2} + (2 {\_}_{3})x_{1}^{4}x_{3}^{2} \text{ is SOS} \]

The model is ready to be optimized by the solver:

JuMP.optimize!(model)

We verify that the solver has found a feasible solution:

JuMP.primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

We can now obtain this feasible solution with:

value(V)

\[ 2.9958453255028186x_{3}^{2} + 12.286107321105156x_{2}^{2} + 15.718362431619155x_{1}^{2} \]

This page was generated using Literate.jl.